Technical  Report  423 


^ ■ 
CQ  - 

rh  • 
CO 

o 

< 

£ 


LEAST  SQUARES,  ADAPTIVE  LATTICE 

ALGORITHMS 


J.  D.  Pack 
E.  H.  Satorius 

April  1979 


D D C 

•?CT?m  nr?i 

SEP  5 1979 

OKSEumd 

D 


Approved  for  public  release;  distribution  unlimited 


NAVAL  OCEAN  SYSTEMS  CENTER 
SAN  DIEGO,  CALIFORNIA  92152 


79  09  5 002  ■{ 


AN  ACTIVITY  OF  THE  NAVAL  MATERIAL  COMMAND 

RR  GAVAZZI , CAPT,  USN  HL  BLOOD 

Commander  Technical  Director 

ADMINISTRATIVE  INFORMATION 

This  report  was  sponsored  by  the  Naval  Ocean  System  Center’s  Independent  Research 
and  Exploratory  Development  Program  (61 152N-ZROOO,  01/ZR014,  08  1 1 . 632-ZR94). 


ACKNOWLEDGEMENTS 

The  authors  wish  to  thank  J.  Treichler  of  Argo  Systems,  M.  Shensa  of  Hydrotronics, 
and  P.  Reeves  and  J.  D.  Smith  of  the  Naval  Ocean  Systems  Center  for  helpful  discussions 
relating  to  this  report. 


Released  by  Under  authority  of 

RH  Hearn,  Head  DA  Kunz,  Head 

Electronics  Division  Fleet  Engineering  Department 


UNCLASSIFIED 

•SECURITY  CLASSIFICATION  of  THIS  PAGE  flWi#n  Dal#  Entorod) 


I 


/ 


L 


[ REPORT  DOCUMENTATION  PAGE 

READ  INSTRUCTIONS 

BEFORE  COMPLETING  FORM 

3 RECIPIENT'S  CATALOG  NUMBER 

4.  TITLE  (ond  Subtitle) 

| LEAST  SQUARES,  ADAPTIVE-LATTICE  ALGORITHMS 

1 CONTRACT  OR  GRANT  NUMBERf«J 

S PEAFOAMING  ORGANIZATION  NAME  AND  ADDRESS 

Naval  Ocean  Systems  Center 

San  Diego,  CA  92 1 52 

10.  PROGRAM  ELEMENT.  PROJECT,  TASK 
AREA  * WORK  UNIT  NUMBERS 

61 152N-ZROOO/01 , ZR014 

08  1 1 i632-ZR94 

11.  CONTROLLING  OFFICE  NAME  AND  AODRESS 

Naval  Ocean  Systems  Center  ■ 1 

San  Diego,  CA  92152 

12.  REPORT  DATE 

|Aprit3*79  / 

13.  NUMBER  OF  PAGES  / ; 1 

30  / ' ; , 

14.  MONITORING  AGENCY  name  4 AOORESS(/f  dllloront  from  Controlling  Olllco) 

. i __ 

is.  SECURITY  CLASS,  (ot  this  r«4»rt; 

Unclassified 

i I ] 

IS#.  DECLASSIFICATION/1  DOWN  GRADING 
SCHEDULE 

l«.  DISTRIBUTION  STATEMENT 

(ot  this  Rmport) 

DISTRIBUTION  STATEMENT  A j 

Approved  lor  public  release;  j 
Distribution  Unlimited 

17.  DISTRIBUTION  STATEMENT  fol  tho  obotroct  an  farad  In  Bloc*  20,  11  dllfaranf  frooi  Roport) 

Approved  for  public  release;  distribution  unlimited. 

IS.  SUPPLEMENTARY  NOTES 

IS.  KEY  WORDS  (Contlnuo  on  raaaraa  a/da  II  nacaaaarr  and  Idanff ly  by  block  numborj 

noise  canceling 
equalization 
signal  processing 

20.  ABSTRACT  (Contlnuo  on  ravaraa  a/da  II  nacaaaary  and  Idontlty  by  Slack  nuaiSaO 

Recently,  it  has  been  shown  by  Morf,  Lee  and  others  that  least  squares,  adaptive  algorithms  may  be  imple- 
mented in  a lattice  form.  This  result  is  of  considerable  interest  due  to  the  rapid  convergence  characteristics  of  least 
squares  algorithms  as  well  as  the  important  properties  of  lattice  structures  (such  as  high  insensitivity  to  round-off 
noise).  This  report  provides  an  explicit  derivation  of  the  joint  real  process,  scalar  least  squares  lattice  algorithm. 

Also,  a Fortran  subroutine  listing  of  the  algorithm  is  presented. 

DD  1473  COITION  OF  I NOV  SS  IS  OBSOLETE 

S/N  01 02-  LF  -01 4-4601 


UNCLASSIFIED 


SECURITY  CLASSIFICATION  OF  THIS  PAOC  fBti#n  Data  Bntmod) 


SUMMARY 


Recently,  it  has  been  shown  by  Morf , Lee,  and  others  that  least  squares,  adaptive 
algorithms  may  be  implemented  in  a lattice  form.  This  result  is  of  considerable  interest 
due  to  the  rapid  convergence  characteristics  of  least  squares  algorithms  as  well  as  the 
important  properties  of  lattice  structures  (such  as  high  insensitivity  to  round-off  noise). 
This  report  provides  an  explicit  derivation  of  the  joint  real  process,  scalar  least  squares 
lattice  algorithm.  Also,  a Fortran  subroutine  listing  of  the  algorithm  will  be  presented. 


Accession  For 

NT  IS  GRA&I 
DOC  TAB 
Unannounced 

1 Justif icatior 

i 

■?,v 

"'5st.ribut.ion/ 

r-bilit.t 

! uvail a 

y 

T Codes 

nd/ or 

' 3.  -A 

1 


CONTENTS 


I.  INTRODUCTION  ...  5 

II.  NOTATIONAL  CONVENTIONS  AND  PRELIMINARIES  ...  7 

III.  DERIVATION  OF  THE  ALGORITHM  ...  9 

III.  A Least  Squares,  One-step  Prediction  Lattice  Algorithm  ...  9 
III- B Joint  Process,  Least  Squares  Lattice  Algorithm  ...  18 

IV.  CONCLUSIONS  AND  DISCUSSION  ...  22 
REFERENCES ...  23 

APPENDIX  A:  DERIVATION  OF  THE  TIME  UPDATE  EQUATIONS  IN 
THE  LEAST  SQUARES  LATTICE  ALGORITHM  ...  25 

APPENDIX  B : FORTRAN  SUBROUTINE  LISTING  OF  THE  LEAST  SQUARES 
LATTICE  ALGORITHM  ...  29 


I.  INTRODUCTION 


In  many  applications,  such  as  speech  processing,  spectral  estimation,  noise 
cancellation,  and  data  equalization,  one  is  interested  in  the  design  of  digital  adaptive 
filter  structures  whose  coefficients  are  continuously  estimated  from  the  incoming  data 
in  such  a way  as  to  minimize  the  mean  squared  error  between  the  filter  output  and  some 
desired  output  sequence.  One  particular  adaptive  filter  implementation  which  has  found 
widespread  use  due,  in  part,  to  its  computational  simplicity  is  the  so-called  least  mean 
squares  (LMS)  adaptive  filter  developed  by  Widrow  (Ref.  1).  This  filter,  which  is  in 
tapped  delay  line  form,  employs  a noisy  gradient  estimation  algorithm  to  compute  its 
coefficients.  Unfortunately,  in  certain  instances  the  convergence  rate  of  the  LMS  filter 
coefficients  to  their  steady-state  values  can  be  quite  slow.  This  slow  convergence  rate 
problem  can  arise,  for  instance,  when  the  spectrum  of  the  input  to  the  LMS  filter 
has  a large  dynamic  range.  This  problem  has  spurred  considerable  recent  interest  in  the 
development  of  adaptive  filter  algorithms  that  converge  much  faster  than  the  gradient-type 
estimation  algorithms. 

A relatively  wide  class  of  rapidly  converging  adaptive  filter  algorithms  arises  in 
the  context  of  a classical  least  squares  problem:  at  each  time  interval  find  the  set  of 
adaptive  filter  coef.  icients  that  minimizes  the  accumulation  of  the  squared  errors  between 
the  filter  output  and  a desired  output  up  to  that  time.  The  extremely  rapid  convergence 
properties  of  the  least  squares  adaptive  filters  have  made  them  appear  promising  in  a 
number  of  different  applications  (e.g..  Refs.  2-5).  One  of  the  major  reasons  for  the 
current  growing  interest  in  least  squares  adaptive  algorithms  is  the  recent  work  of  Morf, 
Ljung,  Lee,  and  others  (Refs.  5-9  and  21)  who  have  shown  how  the  computational 
complexity  of  these  algorithms  may  be  drastically  simplified  from  their  conventional 
implementation  (see,  e.g..  Refs.  2,3).  In  fact,  it  has  been  shown  in  Refs.  5-9  and  21 
that  the  number  of  operations  (multiplications  and  additions)  per  update  for  the  least 
squares  filter  algorithms  can  be  made  proportional  to  the  number  of  filter  coefficients. 

This  is  in  contrast  to  the  conventional  implementation  of  these  algorithms  where  the 
number  of  operations  per  update  is  proportional  to  the  square  of  the  number  of  filter 
coefficients. 

Furthermore,  the  computationally  simpler  least  squares  filters  may  either  be 
implemented  in  tapped  delay  line  (Refs.  4,6)  or  lattice  form  (Refs.  5,  8-9,  and  21). 

The  lattice  realizations  of  the  least  squares  filter  algorithms  are  of  particular  interest 
in  this  report  as  they  offer  a number  of  advantages  over  the  tapped  delay  line 
implementations  (Refs.  5 and  10-12).  Specifically,  in  speech  processing  applications, 
the  lattice  implementations  provide  a simple  check  on  the  stability  of  speech  modelling 
filters  (Refs.  10,  13).  In  all-pole  modelling  applications,  recent  work  (Ref.  14) 
suggests  that  lattice  structures  may  prove  useful  in  determining  the  correct  order  of 
the  model.  In  adaptive  noise  canceling  applications,  lattice  filters  provide  the  capability 
of  being  able  to  dynamically  assign  the  number  of  filter  coefficients  that  proves  most 
effective  at  any  particular  instant  of  adaptation  (Ref.  1 1).  Also,  longer  lattice  filters 
may  be  built  up  from  shorter  ones  by  simply  adding  on  more  lattice  stages.  This 
property  should  prove  useful  in  developing  a “coefficient-slice”  LSI  technology  for 
implementing  adaptive  lattice  filters.  Finally,  an  important  property  of  lattice  filters 
in  general  is  their  high  insensitivity  to  round  off  noise  (Ref.  15). 
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Although  Morf,  Lee,  et  al.,  originally  presented  the  least  squares  lattice  algorithms 
in  Refs.  5,  8-9,  and  21,  the  development  in  these  papers  is  somewhat  limited  and  numerous 
errors  are  present  (primarily  in  Refs.  5 and  8).  It  is  the  purpose  of  this  report  to  provide  a 
more  explicit  (and  therefore  somewhat  lengthier)  development  of  the  least  squares  lattice 
algorithms.  In  particular,  we  will  be  concerned  with  the  pre-windowed,  scalar,  joint  real 
process  lattice  form,  which  proves  useful  in  such  applications  as  linear  prediction,  noise 
cancellation,  and  data  equalization.  Extensions  to  the  vector  input  case  are  straightforward 
and  are  discussed  further  in  Refs.  5,  8,  9 and  21.  We  will  also  present  a Fortran  subroutine 
listing  of  the  least  squares  lattice  algorithm. 


II.  NOTATIONAL  CONVENTIONS  AND  PRELIMINARIES 


Although  the  ideas  which  underlie  the  least  squares  lattice  algorithm  are  simple, 
the  derivation  of  the  algorithm  and  the  algorithm  itsejf  are  quite  technical.  Therefore,  in 
order  to  increase  the  readability  of  our  presentation,  we  will  first  introduce  the  following 
notational  conventions: 

(1)  The  upper  case  letters  (F,  A,  B,  C,  D,  O,  Q,  U,  V,  Y)  will  be  used  to  denote 
column  vectors. 

(2)  The  lower  case  Greek  Letters  (7,  a,  p)  as  well  as  the  lower  case  letters  (d,  e,  f, 
k,  q,  r,  v,  w,  x,  y)  will  be  used  to  denote  scalar  quantities. 

(3)  An  upper  case  script  letter  (/?)  and  an  upper  case  Greek  letter  (II)  will  be  used 
to  denote  square  matrices.  Also,  a prime  (')  will  be  used  to  denote  the 
transpose  operation. 

At  this  point,  we  will  discuss  the  basic  problem  of  interest.  Consider  two  data 
sequences  x(t),  y(t)  t=0,l, ...  ,T,  and  define  Yj^ft)  to  be  an  (N+l)-dimensional  vector 
consisting  of  time-delayed  samples  of  y(t),  i.e., 

YN(t)  = (y(t),  y(t-l), ... , y(t-N))'.  (1) 

We  are  interested  in  obtaining  the  least  squares  filter  error  output  sequence  at  time  T,  i.e., 
x(T)  + F'j^(T)  Yj^(T),  where  F^(T)  is  the  (N+l>dimensional  filter  coefficient  vector  that 
minimizes  the  exponentially  weighted  sum  of  squared  errors: 

T 

^ w7'1  [x(t>  + F^T)  YN(t)]2  . 

The  parameter  w is  a real  constant,  0<w<l,  which  allows  the  filter  to  track  slow  time 
non-stationarities  in  the  data.  Typically,  w is  close  to  1.  The  inverse  of  (1-w)  is 
approximately  the  memory  of  the  algorithm*. 

Differentiating  the  above  sum  with  respect  to  the  components  of  Fjyj(T)  and 
setting  the  result  to  zero  leads  to  the  following  equation  for  the  F^(T)  vector: 

T 

Rn(T)  Fn(T)  = - ^ w7'1  x(t)  YN(t)  . (2) 

In  Eq.  (2),  the  (N+l)  X (N+l)  matrix  fljq(T)  is  given  by: 

T 

RN(T)  = ^ wT_t  YN«)  Yn(‘>  ’ (3) 

t=0 


•It  should  be  noted  that  the  introduction  of  the  exponentially  weighted  sum  of  squared  errors  above 
differs  from  the  treatment  used  in  Ref.  5 to  include  a tracking  parameter  in  the  least  squares  lattice 
algorithms. 
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The  solution  of  Eqs.  (2)  and  (3)  provides  the  least  squares  filter  coefficient  vector,  F^IT), 
at  the  data  sample. 

Two  points  are  worth  noting  concerning  the  above  equations.  First  these  equations 
may  be  applied  to  a number  of  situations.  For  instance,  in  linear  prediction  applications 
y(t)  is  a delayed  version  of  x(t).  In  noise  canceling  applications,  x(t)  is  the  primary  input 
sequence  containing  both  signal  and  noise.  The  sequence  y(t)  is  a noise  reference  input 
that  is  used  to  cancel  the  noise  from  x(t).  In  channel  equalization  applications,  y(t)  repre- 
sents the  received  data  from  the  channel,  and  x(t)  is  a reference  sequence  used  to  train  the 
algorithm.  A second  point  concerning  the  above  equations  is  that  the  limits  on  all  the 
summation  signs  extend  from  t=0  to  t=T.  The  lower  limit,  therefore,  imposes  the  assump- 
tion on  the  data  that  y(t)=0  for  t— 1, ...  , -N.  These  limits  lead  to  the  so-called 
“pre-windowed”  least  squares  algorithm.  If  the  limits  were  t=N  and  t=T,  the  un-windowed 
or  “covariance”  algorithm  is  obtained,  and  if  the  limits  were  t=0  and  t=T+N  then  the  “pre-” 
and  “post-”  windowed  algorithms  would  be  obtained.  A more  complete  discussion  of  the 
different  windowing  methods  may  be  found  in  Ref.  9,  where  a least  squares,  one-step 
predictor  covariance  lattice  algorithm  is  presented. 


III.  DERIVATION  OF  THE  ALGORITHM 

In  order  to  derive  the  least  squares  lattice  algorithm  that  solves  Eq.  (2)  and,  therefore, 
generates  the  sequence  x(T)  + Fj^fT)  Yjg(T)  we  will  start  in  Section  III.  A by  considering  a 
least  squares  prediction  problem.  Then  in  Section  III.B  we  will  show  how  the  solution  of 
the  prediction  pritiem  also  provides  a solution  of  the  basic  problem  of  interest,  i.e.,  Eq.  (2). 

III. A LEAST  SQUARES,  ONE-STEP  PREDICTION  LATTICE  ALGORITHM 

In  deriving  the  least  squares,  one-step  prediction  lattice  algorithm,  we  will  make 
considerable  use  of  some  basic  properties  of  the  Kn(T)  matrix  (n=0,l  ,...,N),  which  may 
easily  be  derived  from  its  definition  in  Eq.  (3)  (with  N replaced  by  n).  Specifically,  we  have 
that: 


Rn(T)  = 


qn(T) 

Q^(T) 

Qn(T)  ; 

/?n-l(T-Dj 

*n-l<T> 

1 Vn(T) 

V„(T) 

! VT>  J 

Rn(  T-D 

+ Y„(T)  Y' l 

as  well  as 


The  dashed  lines  in  Eqs.  (4a)  and  (4b)  denote  matrix  partitioning.  In  these  equations,  the 
scalars  qn(T)  and  vn(T)  are  given  by: 


qn(T)  = wT_t  y2(t)  , 


I 

vn(T)  = ^ WT<  y2(t-n)  , 

and  the  n-dimensional  vectors  Qn(T)  and  Vn(T)  are  given  by: 
T 

Qn(T)=  £ wT-ty(t)Yn.,(T-l), 


V)  - X 


wT-1  y(t-n)  Y n_, (t) . 


In  this  section  we  will  be  interested  in  the  following  least  squares  problem:  find  the 
n -dimensional  coefficient  vector  An(T)  that  minimizes  the  exponentially  weighted  sum  of 
squared  errors. 


T 

Y wT-‘ ej-(t.T)  , 

t=0 

where  en(t,T)  is  the  n1*1  order  “forward”  prediction  evor  residual: 

en(t,T)  = y(t)  + An(T)Yn.1(t-l).  (7) 

The  coefficient  vector  that  minimizes  the  above  sum  is  termed  the  one-step,  (exponentially 
weighted)  least  squares  “forward”  predictor  of  order  n.  As  will  be  seen  in  Section  1I1.B.  the 
problem  of  obtaining  An(T)  or  en(T,T)  is  intimately  related  to  the  problem  of  solving  Eq.  (2). 

Differentiating  the  above  sum  with  respect  to  the  components  of  An(T)  and  setting 
the  result  to  zero  leads  to  the  following  equation  for  the  An(T)  vector: 


1 

£ wT-'v0.lit-i)v;.|(ni 

t=o 


An(D  = ~Qn(T) 


(8) 


Using  the  definition  of  Rn(T)  [Eq.  (3)1  and  keeping  in  mind  that  Yn.j(-1 ) is  the  zero  vector 
(pre-windowed  case),  we  have: 


T 

£ wT-tYn.1(M)Y'.,(t-l)  = /(n.1(T-l). 

t=0 


(9) 


Therefore,  Eq.  (8)  reduces  to 

/?n-i(T-l)An(T)  + Qn(T)  = °n 


where  On  denotes  the  n-dimensional  zero  vector. 

The  solution  of  Eq.  ( 10)  provides  the  vector  of  one-step,  least  squares  “forward" 
predictor  coefficients,  An(T),  of  order  n.  Substituting  the  solution  for  Afi(T)  back  into  the 
sum  of  the  weighted  squared  prediction  errors  yields  the  following  expression  for  the  minimum 
of  this  sum: 


r„(T)  = min 


^ wT_t  e2(t,' 

ft 


,T) 


1 

= ^ w7-1  y2(t)+ A'n(T)Qn(T).  (U) 

ft 


The  last  equality  in  Eq.  (1 1)  follows  from  Eq.  (10).  Note  that  these  equations  may  be  com- 
bined into  the  single,  augmented  matrix  equation 


*n<T)An<D  = 


'r£(T)\ 

On  , 


(12) 
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where  the  extended,  (n+1  )-dimensional  vector,  An(T),  is  given  by 

An(T)  = / * V (,3) 

0 \An(T) ) 

Equation  (12)  follows  directly  from  Eqs.  (4a),  (5a),  (10),  and  (1 1). 

By  analogy  with  the  nth  order,  one-step,  “forward”  prediction  vector,  An(T),  we  can 
also  define  a “backward”  one-step  prediction  vector,  Bn(T),  which  minimizes  the  sum  of 
exponentially  weighted  squared  errors: 


T 

£ wT  t rn(t,T) » 
t=o 

where  rn(t,T)  is  the  nth  order  “backward”  prediction  error  residual: 

rn(t,T)  = y(t-n)  + B^(T)Yn.1(t).  (14) 

As  will  be  seen  in  Section  III.B,  this  “backward”  prediction  vector  will  play  a central  role  in 
the  lattice  formulation  of  the  solution  to  Eq.  (2).  Differentiating  the  above  sum  with  respect 
to  the  coefficients  of  Bn(T)  and  setting  the  result  to  zero  leads  to  an  equation  that  is  analo- 
gous to  Eq.  (8),  i.e., 

r t 

wT-'Yn.,(l)Y^,(t) 

From  the  definition  of  Rn(T)  in  Eq.  (3),  it  can  be  seen  that  Eq.  (15)  may  be  written  as: 

^n-ld)  Bn(T)  + Vn(T)  = On  . (16) 


Bn(T)  = -Vn<T)  • 


(15) 


By  analogy  with  the  minimum  nth  order  forward  squared  error,  r®(T),  we  can  also 
define  a quantity  rf,(T)  which  is  the  minimum  of  the  weighted  sum  of  backward  squared 
errors,  i.e., 


r[j(T)  = min 


wT’t.r2(t,T) 


T 

= £ w7'1  y2(t-n)  + B^(T)  Vn(T)  . 

ft 


(17) 


The  last  equality  in  Eq.  (17)  follows  from  Eq.  (16).  Equations  ( 1 6)  and  (17)  may  be  com- 
bined into  the  single  augmented  matrix  equation: 


Rn(J)  Bn(T)  = 


where  the  (n+l>dimensional  extended  vector  Bn(T)  is  given  by: 


(18) 
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(19) 


Equation  (17)  follows  directly  from  Eqs.  (4b),  (5b),  ( 1 6),  and  (17). 

Another  auxiliary  vector  crucial  to  the  development  of  the  lattice  solution  of  Eq.  (2) 
is  the  (n+1  )-dimensional  vector  Cn(T)  which  is  given  by  the  solution  to: 

fln(T)  Cn(T)  = Yn(T) . (20) 

At  this  point,  we  will  show  how  to  obtain  order  updates  for  the  three  vectors:  An,  Bn,  and 
Cn.  These  order  updates  are  the  main  ingredient  of  lattice  algorithms.  In  particular,  we  will 
first  show  that: 


/An(T)\  kn(T)  / 

WT>  \ ° / ~ rf/T-l)  (®n' 


0 > 

(T-l) , 


(21a) 


where  the  constant,  kn(T),  is  given  by: 
kn(T)  = (last  row  of  Rn+\( T)) 


(r). 


(21b) 


In  order  to  derive  Eq.  (21),  first  note  that  the  right  hand  side  of  Eq.  (21a)  is  in  the 


form: 


where  Dn+]  is  an  (n+1  )-dimensional  vector.  We  will  now  show  that: 

D„t!  * A„+I<T) . «2> 

thereby  establishing  Eq.  (21).  To  verify  Eq.  (22),  premultiply  the  right-hand  side  of  Eq.  (21a) 
by  /?n+|(T)  to  obtain: 


where: 


A 


kn(T)  = (first  row  of  /?n+j(T)) 


(23b) 
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Equation  (23)  follows  from  Eqs.  (4),  (12),  and  (18).  Now,  from  Eq.  (4a)  it  is  seen  that: 

/wT>+Q;,+)(T)Dn+]> 

yQn+l(T) 


*n+1(T)/  1 \=/qn+l(T)  + Qn+‘(T)D"+A 
n‘  \^n+l  / VWT)  + VT-l>Dn+l/ 


^n1 


(T)  - (kn(T)  / rJ,(T-l ))  k 


O 


n+1 


»m), 


where  the  last  equality  follows  from  Eq.  (23a).  Thus,  it  is  seen  that  Dn+j  satisfies: 
*n('r-1>Dn+l  +Qn+1^T)=  °n+l  > 


(24) 


which  is  identical  to  the  equation  that  is  satisfied  by  An+j(T),  i.e.,  Eq.  (10)  with  n replaced 
by  n+1.  Therefore,  assuming  that  J?n(T)  is  nonsingular,*  Eq.  (22)  and,  hence,  Eq.  (21)  are 
verified. 


Note  that  since, 


*n+l(T)  An+1(T) 


n+1' 


/rnH 
"\  0 


iWl<T>\ 


'n+1  / 

we  also  have,  from  Eq.  (23a),  the  following  order  recursion  for  r®(T): 
rn+l^T>  = rn(T>  -M1)  kn(T>  / *n(T*l)  • 


(25) 


(26) 


In  a manner  analogous  to  the  development  of  Eqs.  (21  )-<26),  we  can  also  derive  the 
following  order  update  recursions  for  Bn(T)  and  r{^(T): 

A 


and; 


E (T)  - ° \ kn(T>/^n(T)\ 

Bn+l(T)~  \bn(T-\))  -re(T)  ^ 0 )’ 


rn+|(T>  = rn(T-D  - kn<T>  / rn<T>  • 


(27) 


(28) 


Note  that  in  the  above  order  update  equations,  two  scalars,  kn(T)  and  kn(T),  appear. 
We  will  now  show  that  these  scalars  are  equal.  In  particular,  consider  the  matrix  product: 


/a;(T)  0 \ /An(T)  o \ 

H \ 0 B^(T-l)/  n+,(T)\  0 Bn(T-l)  / * 


(29) 


Notice  that  n is  a symmetrical  2X2  matrix.  Also,  from  Eqs.  (4),  (12),  (18),  (21b),  and 
(23b)  we  have: 


*In  practice,  the  positive  definiteness  (and,  hence,  nonsu.gularity)  of  /?n(T)  can  be  guaranteed  by  initializing 
/?n(T)  to  a positive  scalar  times  the  identity  matrix. 
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’ A^(T) 


ri  = 


ren(T)  kn(T> 

On  °n 


rn(T)  MT) 


(30) 


0 B^T-ir  'kn(T)  rj^(T-l ) ' kn(T)  r{,(T-l)' 


However,  since  n is  a symmetric  matrix,  we  have  the  important  relation: 
kn(T)  = kn(T)  . 


(31) 


Therefore  the  order  updates  [Eqs.  (21a),  (26)  and  (27)  - (28)1  take  on  the  following  form: 

0 \ 


^n+1 


(T)  = 


' An(T) 
l 0 , 


kn(T) 


§n+l(T) 


/ 0 \ ^ 

' ' rnl 


r^d-D  yyT-i^ 
(T)  /An(T)\ 
(T)  \ 0 / 


and, 


ren+l(T)  = ^(T)-k2(T)/rrn(T-l) 


r^l<T)  = rrn(T-l)-kJ(T)/rn(T) 


(32a) 


(32b) 


(32c) 


(3  2d) 


To  complete  the  development  of  the  order  updates,  we  will  now  show  that  Cn(T) 
obeys  the  following  order  recursion:* 


/Cn.,(T)\  rn(T>  - 
Cn(T)=(  +-f-Bn(T). 

n \ 0 / rJ,(T) 


(33) 


To  verify  Eq.  (33),  we  proceed  as  in  the  case  of  Eq.  (21)  and  premultiply  the  right-hand  side 
of  Eq.  (33)  by  7?n(T)  to  obtain: 


Rn(T) 


rn(T)  - 

+ - — Bn(T) 


rrn(T) 


/Yn-l(T)\  / On  ' 

V « /W<T>, 


(34) 


Equation  (34)  follows  from  Eqs.  (4),  (18),  and  (20).  The  constant  a in  Eq.  (34)  is  given  by 
[from  Eqs.  (4b)  and  (20)] : 

a = Vp(T)  Cn.,(T)  = V^(T)  ^-],(T)  Yn_,(T) . (35) 


However,  from  Eqs.  (14)  and  (16)  we  have: 
a = y(T-n)  - rn(T)  . 


(36) 


*ln  Eq.  (33)  and  throughout  the  rest  of  this  report,  we  will  use  rn(T)  in  place  of  rn(T,T)  [see  Eq.  (14)) . 
Likewise,  we  will  use  en(T)  in  place  of  en(T,T)  [see  Eq.  (7)] . 
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Therefore,  Eq.  (34)  reduces  to: 


Rnm 


- VT>  • 


(37) 


Thus,  the  right-hand  side  of  Eq.  (33)  satisfies  the  equation  for  Cn(T)  [Eq.  (20)] . This 
verifies  the  order  recursion  [Eq.  (33)] . 

Equations  (32)  and  (33)  provide  a complete  set  of  order  recursions  for  An(T),  Bn(T), 
Cn(T),  rf^T),  and  r^(T).  From  these  recursions  we  can  obtain  lattice  order  updates  for  the 
residuals  en(T)  and  rn(T).  Specifically,  premultiplying  both  sides  of  Eqs.  (32a)  and  (32b) 
by  V;+|(T)  yields: 


en+1(T)  = en(T)  - (MT)  / rn(T-'))  - (38a) 

and, 

rn+ 1 (T)  = rn(T-l)  - (kn(T)  / r£(T))  en(T) . (38b) 


Likewise,  premultiplying  both  sides  of  Eq.  (33)  by  Yjj(T)  and  defining 

7n(T)  = YJ,(T)  Cn(T)  = Y^(T)  Z?-1  (T)  Y„(T)  , (39) 

we  get  the  following  order  update  relation  for7n(T): 

7n(T)  = 7n-l(T)  + rg(T)  / . (40) 


Equations  (38a)-(38b)  constitute  the  basic  lattice  recursions  that  generate  the  error 
sequences,  rn(T),  which  play  an  important  role  in  solving  our  basic  problem  of  interest,  i.e., 
Eq.  (2).  A schematic  representation  of  these  recursions  is  presented  in  Figs,  la  and  lb.  To 
make  these  recursions  adaptive  in  time,  it  is  necessary  to  derive  a time  update  equation  for 
kn(T).  Such  a derivation  has  been  carried  out  in  Appendix  A,  where  it  is  shown  that  [Eq. 
(A- 1 2)  ] : 


kn(T)  = wkn(T-l)  + 


en(T)rn(T-l) 

>-7n-l(T-l)  ' 


(41) 


With  Eq.  (41),  it  is  now  possible  to  generate  at  each  instant  T the  least  squares  prediction 
error  residuals  en(T),  and  rn(T).  Specifically,  at  each  instant  T we  generate  the  various 
zeroth-order  variables  as  follows: 


e0(T)  = r0(T)  = y(T) , 

(42a) 

rg(T)  = r[,(T)  = wr£(T-l)  + y2(T), 

(42b) 

7_i  (T-l ) = 0 . 

(42c) 
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Then  we  perform  the  various  order  updates  in  the  following  sequence  (n=0, ...  N-l): 

en(T)  rn(T-l) 


kn(T)-  wkn(T',)+  , _7n  |(T-1)  ’ 


(41) 


en+i(T)  = en(T)  - (kn  (T)  / rJ^T-l ))  rn(T-l)  , 


(38a) 


rn+l(T)  = rn(T1  > - (kn(T) ' en<T>  • 


(38b) 


re+|(T)-^r')-k2(T)/rrn(T-l), 

rrn+l‘T‘  = rrn(TI)-kn(T>/rn<T)- 


(32c) 

(32d) 


and,  from  Eq.  (40): 

'fn(T-|)  = >n-l(T-|)  + rS(T-,>/^T-,>- 


(43) 


With  the  generation  of  the  least  squares  residuals  en(T)  and  rn(T)  (n=l, ...  , N)  we  may 
return  to  Eq.  (2)  and  the  generation  of  the  least  squares  error  sequence,  x(T)  + Fj^(T)  Yj^(T). 
However,  before  doing  so,  we  wish  to  make  some  comments  concerning  the  above  least  squares 
lattice  algorithm. 


First,  it  is  somewhat  remarkable  that  the  least  squares  residuals  en(T)  and  rn(T)  obey 
a set  of  order  recursions  which  are  identical  in  structure  to  the  recursions  for  the  minimum, 
mean  squared  error,  one-step  backward  and  forward  linear  prediction  error  residuals  (see, 
e.g..  Ref.  10).  The  latter  recursions  arise  basically  because  of  the  Toeplitz  structure  of  the 
autocorrelation  matrix  associated  with  y(t).  However,  /?n(T)  does  not  have  the  nice  Toeplitz 
matrix  structure,  ar.i  yet  simple-order  recursions  can  still  be  derived  for  en(T)  and  rn(T).  The 
main  reason  for  this  is  that  even  though  /?n(T)  is  not  Toeplitz,  it  can  be  factored  into  two 
Toeplitz  matrices  (see  Ref.  9)  and,  therefore,  still  possesses  “nice”  enough  properties  [i.e., 

Eq.  (4)] , to  allow  a lattice  formulation.  It  should  be  noted  that  this  basic  concept  of  matrices 
that  are  not  Toeplitz  but  nevertheless  are  “close”  to  being  Toeplitz  in  some  sense  has  been 
further  developed  in  the  work  of  Friedlender,  et  al.  (Ref.  16). 


Another  comment  regarding  the  least  squares  lattice  algorithm  derived  above  concerns 
the  parameter  7n(T).  Note  that  this  parameter  only  enters  into  the  lattice  recursions  through 
the  time  update  equation  [Eq.  (41)]  for  kn(T).  In  particular,  the  factor  [ 1 — -yn_ j (T- 1 ) ] — * 
appears  as  a gain  factor,  determining  the  rate  of  convergence  of  kn(T).  An  important  property 
of  7n(T)  is  that  it  is  bounded  by  0 and  1,  i.e., 


0<7n(T)<  1 • 


(44) 


This  can  be  seen  from  the  matrix  inverse  identity  (see,  eg..  Ref.  17): 

. . -R^fT-l)  Yn(T)  Yf,(T)/?~*(T-l) 

/?:'(T)  = w'1  Rl'fT-D-w'2— r— 1 


1 +W-1  Y^(T)  /?~1(T-1 ) Yn(T) 


(45) 
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Substituting  Eq.  (45)  into  Eq.  (39)  gives: 

(w-1  Y' (T)  1 (T-l ) Yn(T)) 

7n(T)= r2 , (46) 

1 +(w->  Y^D^-'d-l)  Yn(T)) 

from  which  Eq.  (44)  is  easily  verified.  Therefore,  when  7n_|(T-l ) approaches  its  maximum 
value  of  unity,  the  factor  ( l-yn.j(T-l  ))"*  becomes  large,  thereby  amplifying  the  gains  in 
the  update  equation  [Eq.  (41)]  for  kn(T).  It  is  interesting  to  note  that  the  main  difference 
between  the  gradient  lattice  algorithms  developed  in  Refs.  10  and  1 1 and  the  least  squares 
lattice  algorithms  is  the  presence  of  the  gain  factor,  ( 1 ^n-i (T-l  ))“*,  in  the  time-update 
equations.  A more  complete  comparison  between  the  gradient  lattice  and  least  squares 
lattice  algorithms  will  be  presented  in  a forthcoming  report  (Ref.  18). 

III.B  JOINT  PROCESS,  LEAST  SQUARES  LATTICE  ALGORITHM 

In  order  to  develop  a lattice  formulation  of  the  least  squares  problem  expressed  by 
Eq.  (2),  we  first  note  that  Eq.  (2)  may  be  rewritten  in  the  form: 

Rn(T)  Fn(T)  + Q£(T)  = On+,  , (47) 


where 


T 

Qn(T)  = J wT't  x(t)  - (48) 

t=0 

and  0 < n < N.  Note  that  in  replacing  N in  Eq.  (2)  by  n in  Eq.  (47),  we  are  actually  con- 
sidering the  larger  problem  of  generating  all  the  least  squares  error  sequences  x(T)  + Fj^T) 
Yn(T),  n=0, ...  , N.  As  will  be  seen,  in  the  lattice  formulation  of  the  least  squares  problem, 
all  of  these  sequences  are  generated  automatically. 

Substituting  Fn(T)  [from  Eq.  (47)]  back  into  the  sum  of  exponentially  weighted 
squared  error  residuals,  i.e., 

T 

^ wT_t  [e*(t,T)|  2 , 
t=0 

where 

ex(t,T)  = x(t)  + F^(T)  Yn(t)  , (49) 

we  obtain  the  following  expression  for  the  minimum  of  this  sum: 

/ T \ T 

Pn(T)  = min  (2  wT‘*  [e£(t,T)l  2 ] = Y w1*1  x2(t)  + F^(T)  Q£(T)  . (50) 

\ t=0  / t=0 

Equations  (47)  and  (50)  may  be  combined  into  the  single,  augmented  matrix  equation: 
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(52) 


where  /?n(T)  is  an  (n+2)  x(n+2)  dimensional  matrix  given  by 


T 

*n(T)  = 2 wT ^ ^n(0  *n(t)’ 

t=0 


and  Yn(t)  is  an  (n+2)-dimensional  column  vector  given  by: 

Yn(t)  = (x(t),  y(t),  y(t-l), ... , y(t-n))' . 

Also,  in  Eq.  (51)  the  (n+2)-dimensional  extended  vector,  Fn(T)  is  given  by 
1 


(53) 


Fn(D- 


Fn(T) 


(54) 


By  analogy  with  Eq.  (4),  it  can  be  easily  seen  from  Eq.  (52)  that  Kn(T)  possesses 


the  following  properties: 

1 

-O 

3X 

H 

w- 

QS'(T) 

/?n(T>- 

Q^(T)  ! 

Rnm 

V iCT)  | 

Vp(T) 

V5'(T)  ! 

vn(T) 

A third  property  of  )?n(T)  is  the  time  shift  relation: 
Rn(T)  = w*n(T-l)  + Yn(T)  Yn(T)  . 


(55a) 


(55b) 


(55c) 


In  Eqs.  (55a)  and  (55b),  the  scalars  q£(T)  and  v*(T)  are  given  by 
T 

q£(T)  = ^ wT-t  x2(t)  (55d) 

tK> 
and 

T 

v*(T)  = ^ wT‘*  y2(t-n)  = vn(T) . (55e) 

t=0 

The  last  equality  in  Eq.  (55e)  follows  directly  from  the  definition  of  vn(T)  (Eq.  (5b)].  Also, 
in  Eq.  (55b),  the  (n+l)-dimensional  vector  V*(T)  is  given  by 


(550 


T 

V*(T)=  £ w1**  y(t-n)  Yn.,(t)  . 
M) 


Using  Eqs.  (5 1 ) and  (55)  as  well  as  the  development  in  Section  II1.A,  we  can  now 
derive  a lattice  solution  to  Eq.  (2).  The  key  relation  in  this  lattice  formation  which  links 
Eq.  (2)  with  the  results  in  Section  IN.A  is  the  following  order-update  relation  for  Fn(T): 


WT>  = 


kn(T)  / 0 \ 

■WT)  VWT>/ 


(56) 


where  the  scalar  k*(T)  is  given  by 


k*(T)  = (last  row  of  /?n+i(T)) 


(57) 


The  derivation  of  Eqs.  (56)  and  (57)  is  analogous  to  the  derivation  of  Eq.  (21 ) and  follows 
by  premultiplying  the  right-hand  side  of  Eq.  (56)  by  /?n+j(T)  and  then  using  Eqs.  (18),  (51 ), 
and  (55).  Premultiplying  both  sides  of  Eq.  (56)  by  Yn+|(T)  leads  to  the  following  lattice 
recursion 


en+l(T)  = e*(T)  - (k£(T)  / r^CT))  rn+,(T) , 


(58) 


where  e*(T)in  Eq.  (58)  is  used  to  denote  e*(T,T)  [see  Eq.  (49) ) . In  direct  analogy  with 
Eq.  (41),  a time  update  equation  for  k*(T)  can  also  be  derived.  Such  a derivation  has  been 
carried  out  in  Appendix  A,  and  the  resulting  update  equation  for  k*(T)  is  given  by  [Eq. 
(A-22)J: 


k*(T)  = wk*(T-l)  + 


e£(T)rn+,(T) 
1 -7n(T) 


(59) 


Equations  (58)  and  (59)  together  with  Eqs.  (32c),  (32d),  (38a),  (38b),  (40),  and 
(41 ) represent  the  complete  lattice  algorithm  that  generates  the  desired  least  squares  error 
sequences  e*(T),  n=0,  ...  , N.  A schematic  representation  of  the  lattice  is  given  in  Fig.  2, 
and  a Fortran  subroutine  listing  of  this  algorithm  is  presented  in  Appendix  B.  It  should  be 
noted  that  in  addition  to  the  zeroth-order  variables  in  Eqs.  (42a)-(42c),  which  must  be 
computed  every  sample  instant  T,  we  also  have 

e^j(T)  = x(T).  (60) 

It  is  interesting  to  observe  that,  as  in  the  case  of  kfi(T),  the  parameter  7n(T)  enters 
into  the  time-update  equation  for  k*(T).  It  is  the  presence  of  this  gain  parameter  which 
enables  the  least  squares  algorithm  to  converge  rapidly.  The  fast  convergence  properties 
of  the  least  squares  lattice  will  be  examined  in  more  detail  in  another  report  (Ref.  18). 
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IV.  CONCLUSIONS  AND  DISCUSSION 


In  this  report,  a complete  derivation  of  the  pre-windowed,  scalar,  joint  real  process 
least  squares  lattice  algorithm  has  been  presented.  This  algorithm,  which  was  originally 
developed  by  Morf,  Lee,  and  others  in  Refs.  5 and  8-9  has  a number  of  useful  properties 
including  quick  convergence  and  a computational  complexity  per  update  that  grows  only 
linearly  with  the  number  of  filter  coefficients.  It  has  been  the  purpose  of  this  report  to 
provide  an  explicit  development  of  the  algorithm. 

Some  comments  regarding  the  application  of  the  least  squares  lattice  algorithm  to 
adaptive  noise  filtering  and  data  equalization  are  worth  making.  First,  note  that  in  the 
present  formulation  of  the  algorithm,  the  residual,  eft(T)  = x(T)  + Fjsj(T)  YN(T),  is 
computed  at  each  iteration.  This  is  in  contrast  to  usual  adaptive  noise  filtering  and 
equalization  algorithms  (see,  e.g..  Refs.  4 and  19),  where  the  residual,  eft(T)  = x(T)  + 
F'N(T-1 ) Yn(T),  is  computed  every  sampling  instant  T.  It  is  straightforward  to  derive  an 
alternate  form  of  the  least  squares  lattice  algorithm  that  computes  eft(T)  instead  of 
efsj(T).  This  alternate  lattice  form  can  be  more  readily  implemented  for  purposes  of 
decision-directed,  adaptive  data  equalization,  as  will  be  discussed  further  in  Ref.  18. 

Another  comment  regarding  the  lattice  algorithm  considered  in  this  report 
concerns  the  choice  of  the  data  window  used  in  developing  this  algorithm.  In  particular, 
we  have  considered  a growing-fading  window,  i.e.,  we  have  exponentially  weighted  the 
data  between  t=0  and  t=T.  The  fade  factor,  w,  allows  the  filter  to  adapt  to  nonstationar- 
ities  in  the  data  as  discussed  previously.  Another  common  method  of  allowing  the  filter 
to  track  data  non-stationarities  is  to  use  a fixed-nonfading  memory  (Ref.  10),  i.e., 
to  uniformly  weight  the  data  between  t=T-TQ+l  and  t=T.  The  parameter  TQ  represents 
the  memory  size  for  the  fixed-nonfading  memory  method.  This  latter  method  may  prove 
useful  for  purposes  of  adaptive  noise  filtering  or  data  equalization  in  highly  nonstationary 
environments  and  is  an  important  subject  for  future  investigation.  See  Ref.  20  for  an 
application  of  this  method  to  intrusion-detection. 
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APPENDIX  A 


DERIVATION  OF  THE  TIME  UPDATE  EQUATIONS  IN  THE 
LEAST  SQUARES  LATTICE  ALGORITHM 

In  this  appendix,  we  will  provide  derivations  of  the  time-update  equations  that 
appear  in  the  least  squares  lattice  algorithms.  In  particular,  we  will  derive  the  time-update 
equations  for  kn(T)  and  k*(T)  in  Sections  A.  1 and  A.2,  respectively. 


A.l  TIME  UPDATE  FOR  kn(T) 

From  Eqs.  (4b)  and  (21b),  we  have  that: 

kn(T)  = v;+i(T)  An(T) . (A- 1 ) 

To  derive  a time  update  equation  for  kn(T),  we  will  separately  derive  time  update  equations 
for  Vj^jd)  and  An(T).  The  update  equation  for  V^+1(T)  follows  directly  from  Eq.  (4c) 
by  noting  that  the  last  row  of  Rn+j(T),  i.e.,  Vn+1(T),  must  obey  the  same  time-update 
relation  as  the  rest  of  the  elements  of  /?n+j(T).  Therefore,  from  Eq.  (4c)  we  see  that 

Vn+i(T)  = wV;+1(T-l)  + y(T-n-l)  Y^T) . (A-2) 

To  derive  a time  update  equation  for  An(T),  we  will  first  consider  An(T).  Note  from 
Eq.  (10)  that 

An(T-l)  = -/?-!1(T-2)Qn(T-l).  (A-3) 


Comparing  Eqs.  (4a)  and  (4b),  it  is  seen  that  Qn(T)  obeys  a time-update  equation  similar  to 
Eq.  (A-2),  i.e.,  from  Eqs.  (4a)  and  (4c): 


Qn(T)  = wQn(T-I)  + y(T)  Yn.,(T-l)  . 

Also,  using  the  matrix  inverse  identity  [see  Eq.  (45)) , it  is  seen  that 

/?n!j  (T-l ) Yn- 1 (T- 1 > Yn- 1 (T-l ) V7'1  > 


(A-4) 


R“!,(T-2)  = w 


K“!l(T-l)  + 


= w 


*n!,(T-l)  + 


1 ~ Yn-,(T-1 ) /?~|,(T-1 ) Yn_j(T-l) 
Vl^-')  Y^., (T-l )/?"!, (T-l) 


l-7n.,(T-l) 


(A-5) 


The  last  equality  in  Eq.  (A-5)  follows  from  Eq.  (39).  Substituting  Eqs.  (A-4)  and  (A-5)  into 
F,q.  (A-3),  we  obtain  (after  combining  terms): 


RnVT-')Yn-l^T-n 

VT-D-VT)*  , -7n.,(T-l) 

- A CTt  + C"-l(T-1)e"(T) 

" l-7n-l(T’1>  ' 


y(T)  + Y^.,(T-1)  An(T) 
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(A-6) 


The  last  equality  in  Eq.  (A-6)  follows  from  the  definition  of  en(T)  [i.e.,  en(T,T)  in  Eq.  (7)] 
as  well  as  Eq.  (20).  Therefore,  from  Eq.  (A-6)  we  have  the  following  update  for  An(T): 

-n(T)  _ . 

(A-7) 


An(T)  = A 


-n' l)  / 0 \ 

" '-WT-ll(cn.,(T.|)J  ’ 


Substituting  Eqs.  (A-7)  and  (A-2)  into  Eq.  (A-l ) gives 

wc  (T)  /S 

kn(T)  = wkn(T-l)  - f.^fr.1)  Vn+l(T-l ) Cn.,(T-l) 

en(T) 

+ y(T-n-l ) Y^(T)  An(T-l)-  ^ } 7n_,(T-l)  y(T-n-l) . 


(A-8) 


In  Eq.  (A-8),  Vn+|(T-1)  denotes  an  n-dimensional  column  vector  whose  elements  are  simply 
equal  to  the  last  n elements  of  Vn+j(T-l).  Note  from  Eq.  (A-2)  that: 

A , A. 


Vp+1(T-1)  = w-1  [V;+1(T)-y(T-n-l)Y^,(T-l)]  . 


(A-9) 


Furthermore,  from  Eqs.  (4a)  and  (4b)  the  following  simple  relationship  may  be  derived: 

Vn+,(T)  = Vn(T-l)  = -/?n.,(T-l)  Bn(T-l) . (A-10) 

The  last  equality  in  Eq.  (A-10)  follows  from  Eq.  (16).  Substituting  Eq.  (A-10)  into  Eq. 

(A-9)  gives 

V;+,(T-l)  = w-1  [-Bjj(T-l) /?n.,(T-l)  -y(T-n-l)  Y^.,(T-1)]  . (A-l  1) 

A 

Equation  (A-8)  may  now  be  simplified  by  substituting  for  V^jfT-l)  from  Eq.  (A-l  1)  and 
expressing  An(T-l)  in  terms  of  An(T)  [from  Eq.  (A-7)] . The  result  is  (after  combining 
terms): 

en(T) 

kn(T)  = wkn(T-l)+T-7n  -(T1)  [y(T-n-l)  + B^(T-l)/?n.1(T-l)Cn.1(T-l)] 

= wkn(T-l)  + en(T)rn(T-l)/(l  -7n.,(T-l)).  (A-l 2) 

The  last  equality  in  Eq.  (A-l  2)  follows  from  Eq.  (20)  as  well  as  the  definition  of  rn(T- 1 ) 

[i.e.,  rn(T-l , T-l ) in  Eq.  ( 14)] . Equation  (A-l  2)  is  the  desired  time-update  equation  for 
kn(T). 

A .2  TIME  UPDATE  FOR  k£(T) 

The  derivation  of  the  time-update  equation  for  k£(T)  follows  along  lines  analogous 
to  those  of  the  derivation  in  Section  A.l . In  particular,  from  Eqs.  (55b)  and  (57) 


kS(T)  = VS;i(T)Fn(T). 


(A- 13) 


26 


The  time  update  for  k*(T)  follows  from  the  time  updates  for  V£+j(T)  and  Fn(T).  From 
Eq.  (55c) 

V*+,(T)  = w V*+,(T-1)  + y(T-n-l)  Y^(T) . (A-14) 

Also,  from  Eq.  (47), 

Fn(T-l)  = -/l“1(T-l)Q*(T-l).  (A-15) 

The  time  update  for  Q*(T)  may  be  obtained  from  Eq.  (55c).  The  result  is: 

Q£(T)  = w Q*(T-I ) + x(T)  Yn(T) . (A- 16) 

Substituting  Eqs.  (A-16)  and  (A-5)  (with  T replaced  by  T+l  and  n replaced  by  n+1)  into 
Eq.  (A-l  5),  we  obtain  (after  combining  terms): 

^j-l(T)  y (T) 

Fn(T-l)  = Fn(T)+  * ~ tx(T)  + Y^(T)Fn(T)] 


Cn(T)e^(T) 


F"(T)+  l-7n(T) 


(A-l  7) 


The  last  equality  in  Eq.  (A-l  7)  follows  from  (20)  and  the  definition  of  e*(T)  (i.e.,  e*(T,T)  in 
equation  (49)).  Therefore,  from  (A- 17)  we  have  the  following  update  for  Fn(T): 


en(T)  / 0 \ 

Fn<T»  = FnCT-„-  r^-m^Cn(T)j  • 

Substituting  Eqs.  (A-18)  and  (A-14)  into  Eq.  (A- 13)  gives 
w e*(T) 

k*(T)  = wk*(T-l)-  Vp+1(T-I)  Cn(T) 


(A-18) 


y(T-n-l)  e5(T)  _ 

+ y(T-n-l)  Y^(T)  Fn(T-l)  - ] -(T)  - Y^T)  Cn(T) . (A-19) 

A 

In  Eq.  (A-19),  V*+](T-1)  denotes  an  (n+l)-dimensional  column  vector  whose  elements  are 
simply  equal  to- the  last  n+1  elements  of  V£+j(T-l).  From  Eqs.  (55a)  and  (4b)  the  following 
simple  relationship  may  be  derived: 

^n+|(T)  = Vn+I(T)  = -Rnm  Bn+,(T).  (A-20) 

The  second  equality  in  Eq.  (A-20)  follows  from  Eq.  (16).  Also,  from  Eq.  (A-14)  we  have 
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Vn+,(T-1) = W>  [V*+,(T)  -y(T-n-l)  Y^T)] 

= -W1  [B;+1(T)/?n(T)  + y(T-n-l)Y^(T)l  , 


(A-21) 


where  in  the  last  equality  we  have  substituted  from  Eq.  (A-20)  for  V*+](T).  Substituting 
Eqs.  (A-21)  and  (A- 18)  into  Eq.  (A-19)  yields  the  following: 

e*(T) 

k*(T)  = wk£(T-l)  + — — [y(T-n-l)  + B;+1(T)/?n(T)Cn(T)l 
en(T)rn+1(T) 


= wk*(T-l)  + 


1 -7„(T) 


( A-22) 


The  last  equality  in  Eq.  (A-22)  follows  from  Eq.  (20)  and  the  definition  of  rn+|(T)  [i.e., 
rn+i(T,T)  in  Eq.  (14)].  Equation  (A-22)  is  the  desired  time-update  equation  for  k*(T). 


APPENDIX  B 

FORTRAN  SUBROUTINE  LISTING  OF  THE  LEAST  SQUARES  LATTICE  ALGORITHM 


The  following  is  a complete  Fortran  subroutine  listing  of  the  least  squares  lattice 
algorithm  given  by  Eqs.  (41),  (42),  (38a),  (38b),  (32c),  (32d),  (43)  and  (58)  through  (60). 

SUBROUTINE  LSANCtU.IP.NTF > 

C BASED  ON  MORF'S  PRE-UINDOUED  FAST-TRACKING  ALGORITHM 

DIMENSION  E<50).P<50).K<50>.RE<S0>,G(50>,RR<S0),EX<S0),KX<50> 
COMMON  'RDATA  ' X<500), Y(500) 

REAL  K fCX  KX& 

NOTE  THAT'j  TAKES  ON  OALUES  1, 2.3, 4,..  UHIIE  N ASSUMES  0. 1,2,3.... 
THE  FOLLOWING  SHOWS  THE  CORRESPONDENCE  BETWEEN  THE  NAMES  USED  IN 
THIS  PROGRAM  AND  THE  NOTATION  FOUND  IN  THE  DERIVATION. 

RECJ)  IS  RE( J#T) 

RR( J } IS  RR( J,T) 

ECJ5  IS  E(J.T) 

TR  IS  RiJ,T-l> 

TPR  IS  RR(J,T-1) 

GMM  IS  GANNA < J-l ,T-1  ) 

GM0  IS  GANNA (J-l.T) 

GOO  IS  GANNA ( J , T ) 


C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


RG  IS  RCJ+l.T) 

RRQ  IS  RRtJ+l.T) 

WRITE  (6,66) 

66  FORMAT  ('1  N T',5X, 'K',8X, 'AKE',7X, 'AKR' ,8X, 'E' ,9X, 'R',8X, 
C RE ',8X,'RR',SX,' GAMMA', 8X,'XICP',8X, 'EX') 

DEL  IS  A SMALL  NUMBER  USED  TO  PREVENT  DIVISION  BY  ZERO. 

DEL5 .OOO0O01 
DO  7 1*1, IP 
R(I )*0. 

RR(I )*DEL 
KXII i-0. 

K;I)*0. 

7 P.ECl  i*DEL 

k:<0*o. 

DO  100  NT-1, NTF 

NTM-NT-i 

TRR-RE(l) 

PEC1  )*U*REUHY<NT)**2 
RRfl  ) *RE( 1 ) 

El  1 1=7 (NT) 

TP-RU  ) 

Ri  I >»V(NT) 

KX0-XCNT )*YCNT )+KX0*W 

EXl 1 )-X(NT)-KX0*Y<NT>/RRU) 

GMM-0. 

GM0-0. 

DO  200  J-I, IP 

N-J-l 

Jl-JM 

C K(J,T)«UtK( J,T-1)+R(J,T-1)*E(J,T)/(1.-GAMMA(J-1,T-1)) 
K(J)-U«(J)+TR*E(J)/(1.-GMM) 

AKE-Kt  J )/RE( J ) 

AKR*K<  J i /TRR 

C ECJ+l.T )*E( J,T )-K( J,T)*R(J,T-1 )/RR(J,T-l) 

E(J1)-E<J)-K(J)*TP/TRR 

C P( J+l ,T  i = R( J,T-1  )-K(J,TUE(J,T)/RE(J,T) 

RO-TP-Kl J )*E'  J i/PEi J > 

C RE(J+l,T)-'RE(J,T)-K(J,T>X*2/RRCJ,T-l) 

RE(J 1 )*REf  J )-<<  J '#<2'TRR 
C PR(J*i,T)«RR<J,T-l)-KiJ,T>**2/RE<J,T) 

RRQ»TRR-K<J)**2/RE<J) 
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C GAMHA(J.T)«GAmiA<J-l<THR<J,T)«a/RR<J,T> 

G00*GN0+R(J)*R(J)/RR(J) 

C KXCJ,T)-U*KX<J.T-l)+EXU,T)*R<J+l,T)/U.-CANNA<J,T>> 

KX l J ) * UtKX  ( J ) *EX  < J ) *R0/  < 1 . -CM  ) 

C EXCJ+1,T)«EX<J,T)-KX<J,T>*R<J+1,T>/RR(J-H,T) 

EX( J1  )*EX( J )-KX( J)*RQ/RRO 
XKP»KX(J)/RRO 

C SUAP  VARIABLES  IN  PREPARATION  FOR  NEXT  ITERATION  ON  J. 
TR*RC J+l ) 

R( J+l )*RO 
TRR«RR<J*1) 

RR( J+l  )»RRO 
GIW-GUl) 

G(J)*GM0 

GNO-GOO 

200  URITE  t6,6  )N,NTN,1C(  J ),  AICE#AKR,E(  J ),R<  J ),RE(  J ),RR(  J ), 
CGM0,XKP,£X<  J ) 

S FORMAT  (2I5,10F10.6) 

10O  CONTINUE 
RETURN 
END 


